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Abstract 

The Lagrange mesh method is a very simple procedure to accurately solve eigenvalue problems 
starting from a given nonrelativistic or semirelativistic two-body Hamiltonian with local or nonlo- 
cal potential. We show in this work that it can be applied to solve the inverse problem, namely, 
to find the equivalent local potential starting from a particular bound state wave function and 
the corresponding energy. In order to check the method, we apply it to several cases which are 
analytically solvable: the nonrelativistic harmonic oscillator and Coulomb potential, the nonlo- 
cal Yamaguchi potential and the semirelativistic harmonic oscillator. The potential is accurately 
computed in each case. In particular, our procedure deals efficiently with both nonrelativistic and 
semirelativistic kinematics. 
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I. INTRODUCTION 



The Lagrange mesh method is a very accurate and simple procedure to compute eigen- 
values and eigenfunctions of a two-body Schrodinger equation . It is applicable for 
both local and nonlocal interactions and also for a semirelativistic kinetic operator, i.e. 
the spinless Salpeter equation . In this method, the trial eigenstates are developed in 
a basis of well-chosen functions, the Lagrange functions, and the Hamiltonian matrix ele- 
ments are obtained with a Gauss quadrature. Moreover, the Lagrange mesh method can be 
extended to treat very accurately three-body problems, in nuclear or atomic physics [3, Q] • 

In this work, we apply the Lagrange mesh method to solve the inverse problem for bound 
states: starting from a given bound state - wave function and corresponding eigenenergy -, 
we show how to compute the equivalent local potential. To our knowledge, this application 
of Lagrange mesh method has not been studied before. It can then be used to compute 
the equivalent local potential of a given nonlocal potential. The determination of equivalent 
local potentials is of particular interest in nuclear physics (see for example Ref. Qj). The 
more interesting point is that our procedure allows to deal with semirelativistic kinematics. 

Our paper is organized as follows. In Sec. [TTl we recall the main points of the Lagrange 
mesh method and show how to apply it to solve a bound state problem with a central 
potential. Then, we give a procedure to compute the equivalent local potential with this 
method starting from a given spectrum in Sec. Illli In order to check the efficiency of our 
method, we apply it to several cases in which the spectrum is analytically known. Firstly, we 
consider three central potentials with a nonrelativistic kinematics in Sec. IIVI the harmonic 
oscillator (Sec. IIV A[) . the Coulomb potential (Sec. IIVBI) . and the nonlocal Yamaguchi po- 
tential (Sec. IIV Secondly, in Sec. El we consider the case of the semirelativistic harmonic 
oscillator for two massless particles, whose solution is also analytical. The accuracy of the 
method is checked in all those cases, and conclusions are drawn in Sec. IVII 



II. LAGRANGE MESH METHOD 

A Lagrange mesh is formed of mesh points Xj associated with an orthonormal set of 
indefinitely derivable functions fj{x) on an interval [a, b]. A Lagrange function fj{x) vanishes 
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(1) 



The weights Aj are hnked to the mesh points Xj through a Gauss quadrature formula 

/.ft N 

/ g{x)dx ^^\kg{xk), (2) 
fe=i 

which is used to compute all the integrals over the interval [a, b]. 

As in this work we only study radial equations, we consider the interval [0, oo[, leading 
to a Gauss-Laguerre quadrature. The Gauss formula ([2]) is exact when g{x) is a polynomial 
of degree 2N — 1 at most, multiplied by exp(— x). The Lagrange-Laguerre mesh points Xi 
are then given by the zeros of the Laguerre polynomial Ln{x) of degree p,J. An explicit 
form can be derived for the corresponding regularized Lagrange functions 

f,{x) = i-iyx-'^' x{x - x.y'L^ix) e-^^\ (3) 

They clearly satisfy the constraint ([T]), and they are orthonormal, provided the scalar prod- 
ucts are computed with the quadrature ([2]). Moreover, they vanish in x = 0. 

To show how these elements can be applied to a physical problem, let us consider a 
standard Hamiltonian H = T{p'^) + V{r), where T{p'^) is the kinetic term and V{r) a radial 
potential (we work in natural units h = c = 1). The calculations are performed with trial 
states I?/)) given by 

N 

\^) = Y.^k\fk) , (4) 

k=l 

where 

(r|/.) = ^^^W^,¥^). (5) 
y nr 

C. is the orbital angular momentum quantum number and the coefficients Ck are linear 
variational parameters, /i is a scale parameter chosen to adjust the size of the mesh to the 
domain of physical interest. If we define r = hx, with x a dimensionless variable, a relevant 
value of h will be obtained thanks to the relation h = Va/xN, where x^ is the last mesh 
point and is a physical radius for which the asymptotic tail of the wave function is well 
defined. This radius has to be a priori estimated, but various computations show that it has 
not to be known with great accuracy, since the method is not variational in h p, llOj ■ 
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We have now to compute the Hamiltonian matrix elements. Let us begin with the po- 
tential term. Using the properties of the Lagrange functions and the Gauss quadrature 1^, 
the potential matrix for a local potential V{r) is diagonal. Its elements are 

POO 

V^, = / dxfiix)Vihx)f,{x) ^ Vihx,) 6^j, (6) 

and only involve the value of the potential at the mesh points. As the matrix elements are 
computed only approximately, the variational character of the method cannot be guaranteed. 



But the accuracy of the method is preserved The matrix elements for a nonlocal 

potential W{r,r') are given by 4| 



CXD 



Wij = h dx dx'fiix) W{hx, hx') fj{x') ^ h ^f\~\j W{hxi, hxj). (7) 



^0 
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The kinetic energy operator is generally only a function of p . It is shown in Ref. 
that, using the Gauss quadrature and the properties of the Lagrange functions, one obtains 
the corresponding matrix 



(p% = ^ 



Xt 



(8) 



where 



(-1)* ^{xiXj) ^/^{xi + Xj){xi- Xj) 2 (i^j) 



(12xf)-i[4+(4iV + 2)xi-xf] («=j)- 
Now, the kinetic energy matrix T(p^) can be computed with the following method j^: 

1. Diagonalization of the matrix p^. If is the corresponding diagonal matrix, we have 
thus = SD'^S"^, where S is the transformation matrix. 

2. Computation of T{D'^) by taking the function T of all diagonal elements of D^. 

3. Determination of the matrix elements Tij in the Lagrange basis by using the transfor- 
mation matrix S: T{p'^) = ST{D'^) S"^. 

Note that such a calculation is not exact because the number of Lagrange functions is 
finite. However, it has already given good results in the semirelativistic case, when T(p^) = 
a/p ^ + m? or even when T{p'^,r) = \/fP -\r U'^{r) 6]. 

The eigenvalue equation H = E reduces then to a system of mesh equations, 

N 

^ [Tij + Vij - E5ij\ Cj = with Cj = ^/hX'ju{hxj), (10) 



where u{r) is the regularized radial wave function and V the local or nonlocal potential 
matrix. The coefficients Cj provide the values of the radial wave function at mesh points. 
But contrary to some other mesh methods, the wave function is also known everywhere 
thanks to Eq. 

III. BOUND STATE EQUIVALENT LOCAL POTENTIAL 

In the previous section, we applied the Lagrange mesh method to solve the eigenequation 
for two-body central problems. We now show that this method allows to solve very easily 
the inverse problem, that is, starting from particular wave function and energy E, to 
find the corresponding equivalent local potential for a given kinematics T. 

In the case of a local central potential, the mesh equations (fTUj) can be rewritten as 



We see from the above equation that, provided we know the radial wave function and the 
energy of the state, the equivalent local potential can be directly computed at the mesh 
points. Let us note that, since the matrix elements depend on the orbital angular 
momentum i, this quantum number has to be a priori specified. The calculation is done 
easily because the potential matrix for a local potential V{r) is diagonal and only involves 
the value of the potential at the mesh points, as shown in Eq. i^. Obviously, this method 
does not require a given normalization for the wave function. Moreover, it is also applicable 
for semirelativistic kinematics. 

We can remark that Eq. (ITTl) contains term which are proportional to u{hxj)/u{hxi). 
They may be difficult to compute numerically with a great accuracy when hxi is either close 
to zero or very large. In these cases indeed, the regularized wave function tends towards zero. 
It means that the ffist values of the potential and also the last ones could be inaccurate. 
It is worth mentioning that, for radially excited states, a particular mesh point Xk could 
be such that hxk is a zero of the wave function. In this case, V{hxk) cannot be computed. 
Although very improbable, this problem could simply be cured by taking a slightly different 
value of N or h. 

In order to check the validity of our method, we will consider four cases where the 
eigenvalue problem is analytically solvable for a given potential V^. This will enable us 




(11) 
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to compare the numerically computed points V{hxi) with the corresponding exact values 
V^{hxi). The number 6, defined by 

V{hx,) - V^{hx,) 



5 = max 



3 <i < N -3} , (12) 



V^{hx^) 

is a measurement of the accuracy of the numerical computations. The more 6 is close to 
zero, the more the method is accurate. The first and last two mesh points are - arbitrarily 
- not included in the computation of 6, since they can introduce errors which are not due 
to the method itself, but rather to a lack of precision in the numerical computations, as we 
argued previously from inspection of formula (fTTj) . 

IV. NONRELATIVISTIC APPLICATIONS 

The kinetic operator which will be used in all the computations of this section is given 

by 

T(P-) = (13) 

where fi is the reduced mass of the studied two-body system. 

A. Harmonic oscillator 



The spectrum of a spherical harmonic oscillator, whose potential reads 

A2 2 

V^^r, = — . (14) 



is given by (see for example Ref. [12, problem 66]) 



Rn,{r) oc /e-^^'/2^1+^/2^Ar2), = A f,-\2n + i + 3/2). (15) 

It is readily computed from the virial theorem that (r^) = (2n + i + 3/2)/ A. Therefore, 
we suggest the following value for the scale parameter: 



h = (16) 

Xn 

4 l {2n + i + 3/2) 
- A ' ^^^^ 

where the factor 4 ensures that the last mesh point will be located in the asymptotic tail of 

the wave function. 
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In order to make explicit computations, we have to specify the value of our parameters. 
We set fi = 0.70 GeV and A = 0.53 GeV^. These parameters can be used in hadron physics 
to roughly describe a cc meson Q. We choose = 30, and the scale parameter is 
computed by using Eq. ( fTTI) . Once these parameters are fixed, Eqs. ( ITTi) and ( |T5i) allow to 
find the equivalent local potential. The result is plotted and compared to the exact harmonic 
potential (|T^ in Fig. [H where we used the wave function in the 25" state (n = 1, £ = 0). 
The numerical result is clearly close to the exact result, and only 30 mesh points are enough 
to provide a good picture of the potential: we have indeed 6 = 2.1 10~^ %, this number 
being computed with Eq. ( fT2|) . The same conclusion holds if other states than the 25* one 
are used, and 6 is always smaller than 1 %. 

In Fig. [21 we show the variation of 6 with the scale parameter h for three different states 
and = 30. We can conclude from this figure that a rather large interval exists where 
the quantity 6 is lower than 1 %. Consequently, the scale parameter does not need to 
be computed with great accuracy: our criterion ( fT6l) is clearly accurate enough since the 
predicted value of h is always located in this interval. The global behavior of 6 which can be 
observed in Fig. [2] is due the difficulty of computing V{hxi) when the scale parameter is too 
small or too large. In this case indeed, the mesh points hxi cover no longer the main part 
of the wave function, and a partial knowledge of the wave function leads to an inaccurate 
description of the potential. 



B. Coulomb potential 

This case is of interest since it enables us to check whether the method we present can 
correctly reproduce a singular potential or not. The radial wave function and eigenenergies 
of a central Coulomb potential 

V"" = -- (18) 



respectively read (see for example Ref. 12, problem 67]) 



R{r) oc r^e--^LltU{2^r), E^^ = (19) 

with Up > 1, < i < Up — 1, and 7 = fin/up. The principal quantum number Up is defined 

by Up = n + i + 1. 
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It can be computed that 




[5nl + 1 - 3e{i + 1)] . 



(20) 



As the evaluation of the scale parameter given by Eq. (1161) yields good results in the harmonic 
oscillator case, it can be adapted to the Coulomb potential, and h is now defined as 



A factor 15 is now needed because the Coulomb potential is a long- ranged one. The wave 
function has thus to be known on a larger domain than for the harmonic oscillator, since 
the latter potential is a confining one. 

In order to numerically compute the equivalent potential from the wave function (ITQll . 
we set yU = 0.70 GeV and k, = 0.27. The particular value of k we chose is commonly used 
in hadron physics to parameterize the one-gluon-exchange part of the potential between 
two heavy quarks [l^. We choose = 30, and the scale parameter is computed by using 
Eq. (EH). The result is plotted and compared to the exact Coulomb potential (ITSl) in Fig. [3] 
for the wave function in the ground state {n = i = 0). The numerical result is close to 
the exact result, with a value of 6 which is equal to 1.4 10~^ %. In particular, the singular 
behavior is well reproduced. To stress this point, we performed another calculation with 
A^ = 100, and h = 0.37 GeV^^ following Eq. ([21]). It can be seen in Fig. [3] that the 
Coulomb potential is then very well matched at short distances. In this case however, we 
have 6 = 0.7 %. Although this precision is still very satisfactory, it seems strange at first 
sight that 6 is higher for a larger number of mesh points. This is due to the fact that the 
mesh points are the zeros of the Laguerre polynomial of degree A^. The first physical point 
which is taken into account in the definition of 6 is hx^, which is smaller for A^ = 100 
{hx3 = 0.811 GeV-i) than for A^ = 30 {hx^ = 0.068 GeV-^). This causes 6 to be larger, 
since the more a point is close to zero, the more the accuracy decreases. 

For what concerns the variation of 6 versus h, the same qualitative features than for the 
harmonic oscillator are observed. Equation ( |2T]) thus appears to give a good evaluation of 
the scale parameter. It can be also checked that a factor smaller than 15 in Eq. ([21]) can 
lead to values of the scale parameter for which 6 is quite larger than 1 %. 




(21) 
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C. Yamaguchi potential 



The Yamaguchi potential is a separable nonlocal potential, given by 

W{r,r') = -v{r)v{r'), (22) 

with 



v{r) = ./p/Jlia + P)e-^\ (23) 

It was introduced in Ref. [15| to study the deuteron (/i = 0.468 GeV). In particular, for 
a = 0.046 GeV and (3 = 0.274 GeV, it admits a bound state whose binding energy is the 
one of the deuteron, that is E = —2.225 MeV. 

A nice particularity of this nonlocal potential is that the bound state wave function can 
be analytically determined. It reads 

-or _ p-0r 

R{r) oc . (24) 

Inserting this wave function into Eq. ffTTl) will provide us with the equivalent local potential 
associated with the Yamaguchi potential. Finding equivalent local potentials coming from 
nonlocal potentials is of interest in nuclear physics, although most studies are devoted to 
scattering states (see for example Refs. [9|]). The bound state equivalent potential of a 
separable nonlocal potential of the form (122|) is shown in Ref. to be given by 

V^(r) = r dr'v(r')u(r'), (25) 
u{r) Jo 

with u{r) the regularized wave function of the bound state for the nonlocal potential. Re- 
lations ( l23l) and (HM can be injected in this last equation to compute that 

V^(r) = - ^. 26) 

As the radial wave function (IMl) is maximal in r = 0, -R(O) (x {jS — a), we can compute 
the scale parameter by demanding that 

R{hxN)/RiO) = e, (27) 

with e a small number, that we will set equal to 10^'^. Then, assuming that a <^ /5 as it is 
the case for the deuteron, h will approximately be given by 

h « _!!LM^. (28) 
9 



The equivalent local potential V^{r) and the one computed with the Lagrange mesh 
method can be compared in Fig. HI The deuteron parameters are used, together with = 30 
and h given by Eq. (1281) . The agreement is satisfactory since 6 = 0.31 %. The extension 
of the wave function is large because the deuteron is weakly bound. An estimation of its 
radius is indeed given by 1.96 fm in Ref. 17|, that is the rather large value of 9.9 GeV~^. 



V. THE SEMIRELATIVISTIC HARMONIC OSCILLATOR 



A nice feature of the Lagrange mesh method is that it allows to solve semirelativistic 



Hamiltonians like the spinless Salpeter equation or the relativistic flux tube model [10|, ll8| , 
which are relevant in quark physics. Equation fllip is consequently applicable if the kinetic 
operator is given by 



T{p^) = 2VF + ^- (29) 
In the ultrarelativistic case where m = 0, the spectrum of the Hamiltonian 



(30) 



can be analytically computed in momentum space in terms of the regular Airy function for 
i = 0. In position space, it reads 13] 



R{r) oc - / dp sin(pr) Ai 



n 



1/3 



p + an 



En 



-(4fi)'/'«n, 



20 



table 10 



where a„ < are the zeros of Ai. They can be found for example in Ref. 

Thanks to the particular properties of the Airy function, it can be computed that 

2/3 



(31) 



13]. 



2l| 



ly^'on 

nj 3 ■ 



(32) 



The scale parameter will thus be computed with the relation 

h 



4 /2V/' I a. 



xn / V 3 

in analogy with the similar case of the nonrelativistic harmonic oscillator. 

The comparison between the potential computed with our method and the exact one 



(33) 



(34) 
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is g iven in Fig. The value Q = 0.2 GeV^ is typical for potential models of light quarks 



131]. But, we present our results as dimensionless quantities. The curves are thus universal: 



they do not dependent on Q, which is the only parameter of this Hamiltonian. Although 
still satisfactory, the agreement is not as good as with the nonrelativistic applications. We 
find indeed 6 = 3.1 %. By inspection of Fig. [5l it can be seen that the last points slightly 
differ from the exact curve. These points are related to the value of the wave function in 
its asymptotic tail, as it can be seen from Eq. ffTTj) . It means that finding the equivalent 
potential, especially with a semirelativistic kinematics, needs a good knowledge of the tail, 
which is not often necessary for computation of the energy spectra. 

In our case, the discrepancies for the last points are due to the computation of the wave 
function in the asymptotic regime. It can be checked that a resolution of Hamiltonian ( 130|1 
with the Lagrange mesh method leads to a wave function which asymptotically decreases 
faster than the exact wave function, given by Eq. (1311) . Conversely, if one starts from the 
exact wave function, the Lagrange mesh procedure will lead to a potential which does not 
increase enough asymptotically, as we observe in Fig.O Fortunately, only the very last points 
are affected, as it is shown in Fig. El By varying and h, that is to say by varying the 
interval where the potential is computed, one can always correctly reproduce the potential 
in a given region: the more hx^ is large, the larger is the interval where the potential is 
correctly reproduced. Finding the equivalent potential with a spinless Salpeter equation 
seems thus to require a more careful study: several curves have to be computed by varying 
h and in order to understand whether the long range behavior of the potential is physical 
or simply due to a numerical artifact. 



VI. CONCLUSIONS AND OUTLOOK 



In this work, we extended the domain of application of the Lagrange mesh method to a 
particular type of problem: to find the equivalent local potential corresponding to a given 
bound state with a given kinematics. We assumed a central problem. Starting from a 
particular radial wave function and the corresponding energy, the method we presented here 
allows to compute the equivalent local potential at the mesh points. We checked the accuracy 
of the computations in various cases whose solutions are analytically known. Firstly, we 
studied the well-known nonrelativistic harmonic oscillator and Coulomb potentials. These 
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potentials are correctly reproduced by the Lagrange mesh method with a precision better 
than 1 %, provided the scale parameter is large enough to take into account the asymptotic 
tail of the wave function. Moreover, the singularity of the Coulomb potential is well matched. 
The numerical parameters are the number of mesh points, and the scale parameter. It 
appears that a typical value of 30 mesh points is enough to provide a good picture of the 
potential. As it was the case for usual eigenvalue problems, the scale parameter does not 
need to be accurately determined: a rather large interval exists where the precision is lower 
than 1 %. 

If the spectrum comes from a nonlocal potential, our method will compute the equivalent 
local potential. This problem is of interest in nuclear physics jol. As an illustration, we 
applied it to the nonlocal Yamaguchi potential describing the deuteron. In this particular 
case, the spectrum is analytical as well as the corresponding equivalent potential. Again, 
the accuracy of our method is very good. 

Finally, our procedure can also be easily adapted to the case of a semirelativistic kine- 
matics. As a check, we studied the semirelativistic harmonic oscillator. Again, the potential 
is correctly reproduced, but it appears that the asymptotic behavior of the potential is prob- 
lematic. This is an artifact of the method in the semirelativistic case: by varying the mesh 
size, one can indeed see that the value of the potential at the last mesh points is systemati- 
cally too low, but the harmonic shape of the potential is well reproduced at the other mesh 
points. 

Our purpose is to apply this method to the study of systems containing quarks and 
gluons. In particular, glueballs, which are bound states of gluons, are very interesting 
systems because their existence is directly related to the nonabelian nature of QCD. Bound 
states of two gluons can be described within the framework of potential models by a spinless 
Salpeter equation with a Cornell potential: a linear confining term plus a Coulomb term 
coming from short-range interactions [22]. Such a phenomenological potential has been 
shown to arise from QCD in the case of a quark- ant iquark bound state 231]. Theoretical 
indications show that it could be valid also for glueballs {zi]. Moreover, recently, the mass 



and the wave function of t 
computed in lattice QCD 



le scalar glueball (with quantum numbers J^'-' = 0~^~^) has been 



25| . Thanks to the Lagrange mesh method, these data could be 



used to extract the potential between two gluons from lattice QCD, and see whether it is a 
Cornell one or not. This study will be published elsewhere. 
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FIG. 1: Comparison between the potential computed from the 25 wave function (llSp thanks to the 
Lagrange mesh formula (jlip . and the exact harmonic potential given by Eq. (jl4p . The equivalent 
potential is only known at the mesh points (circles), and the exact potential is plotted with a solid 
line. The regularized wave function is also plotted with an arbitrary normalization (dashed line). 
We used = 0.70 GeV, A = 0.53 GeV 2, = 30, and /i = 9.8 lO^^ GeV^^ following formula ([H]). 
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FIG. 2: Evolution of 5 versus the scale parameter h for the 15" (full circles), IP (empty circles), and 
25" (triangles) states for N = 30. The gray boxes are the different values of 5 for a scale parameter 
computed with formula ()17p . They all ensure a value of 8 lower than 1 %. 
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FIG. 3: Comparison between the potential computed from the 15 wave function (jl9p thanks to the 
Lagrange mesh formula and the exact Coulomb potential given by Eq. (|18|) . The potential has 
been computed with = 30 (circles) and = 100 (crosses) mesh points, but only a few points are 
plotted for clarity. The regularized wave function is also plotted with an arbitrary normalization 
(dashed line). We used fi = 0.70 GeV and k = 0.27. For N = 30, it can be computed that 
h = 1.3 GeV-^ and 6 = 1.4 10"^ %; for N = 100, we have h = 0.37 GeV~i and 6 = 0.7 %. 
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FIG. 4: Comparison between the equivalent local potential (circles) computed from the wave 
function (j24p with E = —2.225 MeV and the exact equivalent local potential (solid line) given by 
Eq. ()26p . The regularized wave function is also plotted with an arbitrary normalization (dashed 
line). We used a = 0.046 GeV, /3 = 0.274 GeV, and N = 30. Following formula (i28]). h = 
1.8 GeV-i. 
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FIG. 5: Comparison between the potential computed from the IS wave function (I3ip with a 
semirelativistic kinematics (circles) and the exact harmonic potential (solid line) given by Eq. (I34p . 
The regularized wave function is also plotted with an arbitrary normalization (dashed line). We 
used N = 30, and (J^/2)V3/j = 0.034 from Eq. 





FIG. 6: Complete curves giving the potential computed from the IS wave function (j3ip with a 
semirelativistic kinematics for several choices of h and N . These choices ensure more or less the 
same mesh point density. The potential is only known at the mesh points (symbols), and the exact 
harmonic potential is plotted with a solid line. 
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